
dataname<-c("summary","backprob")
eventname<-c("Limitation","Back Problem")
datanum<-1
eventmeans<-TRUE

cutlim<-c(0,0.5)
for(event in c("workprob","backprob")){
    panelnames<- c( "" 
    )
    vars<-c("yetapplyDI_rec", "yetappealDI_rec", "onDI_rec") 
    
    varnames<-c("Applying for DI","Applied for DI", "Reapplied/Appealed", "Receiving DI" 
                )
    
    vargroupnames<-list(c("Applied for DI", "Reapplied/Appealed", "Receiving DI") 
    )

    cutlabs<-list(c("Single","Partnered"),
                  c("Renter","Owner"))
    
    gnum<-1
    for(gvar in c("haspartner",
                  "homeowner")){
      source("./hrs_eventpreamble.R")

      cutlab<-cutlabs[[gnum]]
      decimals<-2
      if(gvar != "sp_wagerank"){
      tab<-  TR(c("","(1)","(2)","(3)","(4)","(5)")) +
        TR(c("",paste0("Pre-",eventname[datanum]," Means"),paste0("Post-",eventname[datanum]," Means")),cspan=c(1,2,3)) + 
        TR(c("","Single","Partnered","Single","Partnered","Partnered - Single"))+
        midrulep(list(c(2,3),c(4,6)))  +
        TR(c("Panel A: Application Behaviors"),cspan=c(6), surround = "\\textit{%s}")
      }
      

      hrsdata[!is.na(group),pval:=feols(as.formula(paste0("as.numeric(as.character(group)) ~ 1 | interaction(time,agebin)")),
                                        data = hrsdata)$fitted.values]
      
      hrsdata[,pweight:=NULL]
      hrsdata[group==1 & pval < 1 & pval > 0,pweight:=1]
      hrsdata[group==0 & pval < 1 & pval > 0,pweight:=pval/(1-pval)]
      hrsdata[time==0&!is.na(pweight),weighted.mean(age,pweight),by=c("group")]
      hrsdata[time==0,mean(age),by=c("group")]
      
      vnum<-1
      for(vvar in vars){
      meandata<-hrsdata[time>= -6& time  <= 20,.(mean(get(vvar)),N=sum(N)),by=c("time","group")]
      meandata<-meandata[order(time),]
      
      
      pdf(paste0("./HRS_eventstudy/",event,"_level_",vvar,"_",gvar,".pdf"))
      print(ggplot(data= meandata[time>=-4 & time <=4 & !is.na(group),])+
              geom_line(aes(x=as.numeric(as.character(time)),y=V1,group=as.factor(group),linetype = as.factor(group)),size=1)
            +geom_vline(xintercept=0,linetype="dashed")
            +geom_hline(yintercept=0,linetype="dashed")
            +ylim(cutlim)
            +scale_x_continuous(breaks=-4:4)
            +scale_linetype_manual(labels=cutlab,values=c("dashed","solid"))
            +theme(legend.position="bottom")
            +labs(x="Years from Onset",y=varnames[vnum],linetype="")
            +theme(legend.key.size = unit(2,"line"))
      )
      dev.off()
      
      
      feols(as.formula(paste0(vvar," ~ as.factor(time)+as.factor(time):group -1")), data = hrsdata[time>= 0 & time <= 4&!is.na(pweight),],cluster="id",
            weights=hrsdata[time>= 0 & time <= 4&!is.na(pweight),]$pweight)
      
      feols(as.formula(paste0(vvar," ~ group")), data = hrsdata[time>= 0 & time <= 4 &!is.na(pweight),],cluster="id",
            weights=hrsdata[time>= 0 & time <= 4&!is.na(pweight),]$pweight)
      means<-feols(as.formula(paste0(vvar," ~ group-1")), data = hrsdata[time>= 0 & time <= 4&!is.na(pweight),],cluster="id",
            weights=hrsdata[time>= 0 & time <= 4&!is.na(pweight),]$pweight)
      premeans<-feols(as.formula(paste0(vvar," ~ group-1")), data = hrsdata[time< 0 & time >= -4&!is.na(pweight),],cluster="id",
                   weights=hrsdata[time< 0 & time >= -4&!is.na(pweight),]$pweight)
      meandata<-hrsdata[time>= -6& time  <= 20 &!is.na(pweight),.(weighted.mean(get(vvar),pweight),N=sum(N)),by=c("time","group")]
      meandata<-meandata[order(time),]
      
      
      pdf(paste0("./HRS_eventstudy/",event,"_agebal_level_",vvar,"_",gvar,".pdf"))
      print(ggplot(data= meandata[time>=-4 & time <=4,])+
              geom_line(aes(x=as.numeric(as.character(time)),y=V1,group=as.factor(group),linetype = as.factor(group)),size=1)
            +geom_vline(xintercept=0,linetype="dashed")
            +geom_hline(yintercept=0,linetype="dashed")
            +ylim(cutlim)
            +scale_x_continuous(breaks=-4:4)
            +scale_linetype_manual(labels=cutlab,values=c("dashed","solid"))
            +theme(legend.position="bottom")
            +labs(x="Years from Onset",y=varnames[vnum],linetype="")
            +theme(legend.key.size = unit(2,"line"))
      )
      dev.off()
      


#####################################################################
#####################PLOTTING RESULTS################################
      
    if(gvar !="sp_wagerank"){
      namepos <- vnum


      tab <- tab + 
        TR(vargroupnames[[1]][namepos])%:%
        TR(c(premeans$coefficients["group0"] , 
             premeans$coefficients["group1"],
             means$coefficients[["group0"]],
             means$coefficients[["group1"]],
             means$coefficients[["group1"]] - means$coefficients[["group0"]]
        ),
        pvalues=c(ifelse(is.na(premeans$coeftable["group0","Pr(>|t|)"]),1,premeans$coeftable["group0","Pr(>|t|)"]) ,
                  ifelse(is.na(premeans$coeftable["group1","Pr(>|t|)"]),1,premeans$coeftable["group1","Pr(>|t|)"]),
                  ifelse(is.na(means$coeftable["group0","Pr(>|t|)"]),1,means$coeftable["group0","Pr(>|t|)"]),
                  ifelse(is.na(means$coeftable["group1","Pr(>|t|)"]),1,means$coeftable["group1","Pr(>|t|)"]),
                  2*pt(abs(means$coefficients[["group1"]] - means$coefficients[["group0"]])/(
                    (means$se[["group1"]]^2+means$se[["group0"]]^2)^0.5),
                    df = means$nobs -means$nparams,
                    lower.tail=FALSE)
        ), dec=decimals
        )+
        TR("")%:%
        TR(c(premeans$se["group0"] , premeans$se["group1"],
             means$se[["group0"]],
             means$se[["group1"]],
             (means$se[["group1"]]^2+means$se[["group0"]]^2)^0.5), se=TRUE)
      vspace(6) 
      
      }
      vnum<-vnum+1
      }

 var <- "wait_time"

  means<-feols(as.formula(paste0(var, "~ group -1")),
               data = hrsdata[time == 0 & !is.na(pweight),],
               weights=hrsdata[time == 0 & !is.na(pweight),pweight], 
               cluster="id")
  if(gvar !="sp_wagerank"){
    tab <- tab + 
      vspace(5)+
      TR(c("Panel B: Statistics on Eventual DI Beneficiaries"),cspan=c(6), surround = "\\textit{%s}")+
      TR(c("Days to secure DI allowance"))%:%
      TR("")%:%TR("")%:%TR(c(means$coefficients["group0"], 
                             means$coefficients["group1"],
                             means$coefficients["group1"] - means$coefficients["group0"]),
                           pvalues = c(means$coeftable["group0","Pr(>|t|)"],
                                       means$coeftable["group1","Pr(>|t|)"],
                                       2*pt(abs(means$coefficients[["group1"]] - means$coefficients[["group0"]])/(
                                         (means$se[["group1"]]^2+means$se[["group0"]]^2)^0.5),
                                         df = means$nobs -means$nparams,
                                         lower.tail=FALSE)), dec = 1
      )+
      TR("")%:%TR("")%:%TR("")%:%TR(c(means$se["group0"], 
                                      means$se["group1"],
                                      (means$se["group1"]^2 + means$se["group0"]^2)^0.5), se = TRUE, dec = 1)
    
  }
  
  TS(tab, file=paste0(paste(c(dataname[datanum],gvar,"pooled_tablemeans"),collapse="_")),
     output_path="./HRS_eventstudy/",
     pretty_rules=T,
     header=c('r',rep('c',5)))
  gnum<-gnum+1
    }
  datanum<-datanum+1
}




